Demographic noise and resilience in a semi-arid ecosystem model 
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^v^j "The scarcity of water characterizing drylands forces vegetation to adopt appropriate survival strategies. Some of these 
.generate water-vegetation feedback mechanisms that can lead to spatial self-organisation of vegetation. To date these 
^^phenomena have mostly been studied with models representing plants by a density of biomass, varying continuously 
jj^ ,in time and space. Such models disregard the discrete nature of plant individuals and their intrinsically stochastic 
behaviour. These features give rise to demographic noise, which can influence the qualitative dynamics of ecosystem 
models. In the present work we explore the effects of demographic noise on the resilience of a model semi-arid ecosystem. 
We introduce a spatial stochastic eco-hydrological hybrid model in which plants are modelled as discrete entities subject 
r ~" l to stochastic dynamical rules, while the dynamics of surface and soil water are described by continuous variables. The 
' model has a deterministic approximation very similar to previous continuous models of arid and semi-arid ecosystems. 
I By means of numerical simulations we show that demographic noise can have important effects on the extinction and 
_0 recovery dynamics of the system. In particular we find that the stochastic model escapes extinction under a wide range 
O of conditions for which the corresponding deterministic approximation predicts absorption into desert states. 
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f— i "1. Introduction 

In arid and semi-arid ecosystems water constitutes the 
H -main limiting resource for vegetation. The harsh envi- 
"ronmental conditions, related to frequent droughts, limit 
OO plant recruitment and survival and often prevent veg- 



OO 

in 



etation from fully covering the ground. In such ar- 
eas vegetation largely occurs in patches of high den- 
sity surrounded by bare soil, and fo r ming^ spatial vegeta- 

o 

(N 



tion p atterns (e.g. lAguiar and Salal . Il999t iBarbier et al 



20061: iDeblauwe et all l2008h . These spatial structures 
emerge from the system's dynamics as a consequence 
of scale-dependent water-vegetation feedback mechanisms 
even in absence of any underlying spatial heterqgene- 



^ iRietkerk et al 



Klausmcicr, 1999; von Hardcnbcrg et al., 2001 



2002tlMeronLl201ll) . For example, in arid 



areas the vegetation increases the infiltration capacity of 
the soil as compared to bare ground because of root pene- 



tratio n (jRietkerk and van de Koppell . 1199 7; Walke r et al 



198lh . While vegetation patches compete for water, vege- 



tation enhances its own growth locally within the patches. 
This so-called infiltration feedback is known to be one 
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of the most important scale-dependent water-vegetation 
feedback mechanisms, leading to self-organization and pat- 



tern formation phenomena (e.g. IRietkerk et all [2002) 



Scale-dependent resource concentration mechanisms of 
this type are also connected to the possible occurrence 
of catastrophic shifts in ecosystems. For example, a veg- 
etated patchy state may turn into a degraded state with 
mostly bare soil if rainfall decreases below a certain thresh- 
old. A subsequent increase in rainfall above the threshold 
may not be enough to recov e r the previous vegetation state 
(e.g. IRietkerk et all 120041: IRietkerk and" van de Koppell 
l2008t iBaudena and Rietkerkl 12012ft . This is an example 
of how the resilience of such ecosystems is strongly inter- 
related with spatial structure. By 'resilience' we mean here 
the ability of an ecosystem to remain in a given domain 
of attraction and to return quickly to the same st ate after 
disturbances ( Rietkerk and van de Koppeil2008h . 



In models used to represent dryland ecosystems both 
water and plants are often represented as density fields, 
varying continuously in time and space. In these models 
their dynamics is typic ally described by deterministic dif- 
feren tial equations (e.g. iGilad et all 120041 ; IRietkerk et al 



2002h . In other models the plant dyna mics is modelled us 



ing s tochastic differential equations 

((m 

anor and Shnerb. 

2008), or instead the vegetation is represen ted by dis- 



crete s tates in a drastically simplified way (e.g. iKefi et al 
2007alibh . 



In drylands, vegetation is composed of a relatively 
small number of plants coexisting with regions of empty 
land in which the density of biomass is essentially zero. 
Under these conditions the intrinsic stochastic behaviour 
of individual plants can b e relevant for under s tandi ng the 
behaviour of the system ([Black and McKand . l2012i ) . The 
spatial self-organization of these ecosystems has not yet 
been explored with models representing single plants in- 
dividually, even though such individual-based approaches 
are now increasingly common in other areas of ecology 
and biology. Demographic stochasticity, due to the 
discreteness of the underlying population, is well known 
to have important consequences on the dynamics. These 
are lost when modelling the sy stems of interest in terms o f 
continuous densities (see e.g. Black and McKanel l2012h . 
Such deterministic approaches are formally valid only in 
the limit of infinite populations. In finite populations 
one finds systematic deviations from the deterministic 
description; this is due to the intrinsic stochastic dynamics 
of the individuals constituting the system. Generally 
speaking, stochastic effects become more pronounced 
as the number of individu a ls in t he population is re 



duced (iRamaswamv et al" 



Biancalani ct al 



Biancala ni et al... _ 
McKane and Newmanl . 



Roger s et al 



2012 



Butler and Goldenfeld, 2011 



Butler and Goldenfeldl . 120091 



20051 ) . In plant ecology, forest 



models are a success f ul example of individual- based mod- 
elling dBotkin et all 1 19721: IShugartl. Il984t IPacala et al 



1993, Il996t iMoorcroft etail [2001 



IPerrv and En right 



20061) . The approa c h has been extended to gra s slands 



as well (lJ0rgensenl. |201 It ICoffin and Lauenroth . 1990: 



Peters! 2002 : Rastetter et all " 2003 ). and it has also been 



used as a method to paramet rize mean-field differen tial 
equations from fine-scale data ( Moorcroft et alll200ll ) 

The consideration of plant individuals naturally leads 
to an individual-based representation. In arid areas vege- 
tation is strictly dependent on water availability, and the 
dynamics of water is effectively captured by a determin- 
istic continuous approach. Thus, a natural approach to 
represent the dynamics of drylands is in terms of so-called 
'hybrid models', in which the dynamics of water is repre- 
sented by differential equations, and in which the plant dy- 
namics are descr i bed b y an individual-based model (IBM) 
( Vincenot et al. . 201dl ). Thus, the model combines both 



continuous and discrete variables. 

In our work we use such a hybrid model to study the 
resilience of semi-arid ecosystems against desertification. 
The model includes the water-vegetation infiltration feed- 
back, and thus we expect vegetation patterns to emerge 
as previously observe d in other type s of m odels of semi- 



arid ecosystems (e.g. iRietkerk et all 120021) . The use of 



individual-based models can be computationally demand- 
ing if too many details are included, and this can ob- 
struct the identification of the underlying mechanisms rel- 
evant for the behaviour of the system. Our approach 
constitutes a compromise between realism, tractability, 
and insight, concentrating on the most relevant vegeta- 



tio n features. We extend an IBM, recently investigated 
by iRealpe-Gomez et al. ( 2012 ) to include spatial interac- 
tions. We study the resilience and stability of the model 
dryland ecosystem when the individual nature of plants 
and its intrinsic stochasticity are considered. Our aim is 
to investigate the role of demographic noise and to under- 
stand whether noise enhances recovery or whether it may 
drive the system to extinction. To do so, we will compare 
the outcomes of a stochastic IBM and a deterministic con- 
tinuous model. 

2. Model definitions 

2.1. Stochastic model 

2.1.1. General considerations 

In this section we introduce our stochastic hybrid model 
for semi-arid ecosystems. The model describes the dy- 
namics of vegetation, soil and surface water in a given 
area to which the model is applied. The source of sur- 
face water is rainfall. Surface water then infiltrates the 
soil where plants can take it up. Plants are consid- 
ered as discrete entit i es fol lowing a stochastic dynamics 
I Realpe-Gomez et all 120121 ). while soil and surface wa- 



ter are modelled by co ntinuous variables whose dynamics 
would be deterministic Jpueyo et all 120081: [Rietkerk et al 



2002; iHilleRisLambers et all 12001) were they not influ- 
enced by plants. Figure [1] summarizes the dynamics of the 
model. For computational simplicity we will work mostly 
in one spatial dimension. Occasionally we will also con- 
sider the case of two dimensions; this is to illustrate that 
the results of the one-dimensional model are still valid in 
this more realistic case. We will define the model in one 
spatial dimension only, the modifications required to ex- 
tend it to two dimensions are straightforward. In the one- 
dimensional model we assume that the land is partitioned 
into L uniform square cells, each labelled by an index i. 
Thus the model describes a linear array of square cells. 
The length of the side of a cell is denoted by h. With ev- 
ery cell i we associate three variables corresponding to the 
discrete (integer) number of plants, rij, in the cell and the 
continuous quantities of soil and surface water, Wj and Uj, 
respectively, on that piece of land. We will assume that 
the density of biomass per unit area in cell i is given by 
Pi = firii, where \x represents the mass of one plant indi- 
vidual divided by the area of the cell, i.e. fi = mp/h 2 . We 
have made here the simplifying assumption that all plants 
have identical mass mp. Thus, \x characterizes the contri- 
bution of an individual plant to the biomass density. We 
consider the cells as homogeneous, i.e. we do not take into 
account any structure within the cell, such as the spatial 
distribution of plants within a cell. We also neglect prop- 
erties of individual plants, such as for example different 
plant sizes or stages of growth. The detailed dynamics of 
the various components of the model are explained below, 
where we also introduce the relevant model parameters. 
The units and numerical values of all the parameters in 
the model are summarised in Table [TJ 
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2.1.2. Water dynamics 

For all cells, i, the dynamics of the depths of soil and 
surface water, u>i and Oj, (measured in mm) are described 
by differential equations that are deterministic when the 
biomass density pi (measured in g mT 2 ) is kept constant; 
namely 



^ = F rJ {p i ,LO i ,(T l ) + D a A(Ti . 



(la) 
(lb) 



In each of these equations the first term on the right-hand 
side describes the water dynamics within a cell. These are 
captured by the functions F^ and F a , to be specified below 
(see Eqs. ([3]) and ([!])). The second term in each of the two 
equations describes spatial transport proc esses; these pro 



cesse s are assumed to be diffusive (see e.g. iRietkerk et al 



2002), and we model them using the discrete diffusion op- 
erator 



Auii 



Act, 



^2 Yl ( w i- w *)» 



J£ Yl fa" *)- 



(2a) 
(2b) 



Here the sum is over all elements j € N(i), i.e. the set of 
cells j which are nearest neighbours of cell i. The coeffi- 
cients D u and D a of the diffusion operator in Eqs. ([1]) are 
the diffusion constants corresponding to the transport of 
soil and surf ace water, respectively. 



Following Hi lleRisLambers et al.l (|2001l) we define the 
functions F u and F a , describing the on-site water dynam- 
ics, as 



F u (p,u},a) = a(p)a - /3{uj) p 
Fo-(p,u),cr) = R - a(p)a . 



rui , 



Here, 



a{p) 



p + k 2 W 



(3a) 
(3b) 



(4) 



p + k 2 ' ^ v ' lj + ki ' 

describe the infiltration rate of surface water into the 
soil, and the soil water uptake due to plants, respectively. 
These rates saturate for large values of p and w, respec- 
tively. The model parameters k\ and k 2 determine the 
exact shape of the saturat ion curves and are known a s 
'half-saturation' constants (jHilleRisLambers et al. . 2001 ). 
The dependence of a on the biomass density, p, introduces 
the infiltration-watcr-vcgctation feedback, known to gen- 
erate spatia l vegetation patterns in existing deterministic 
models fe.g. IRietkerk et all 120021 ). The parameter r char- 
acterizes the loss of soil water due to evaporation (Eq. 
(J3^)), while the parameter R in Eq. ((2p) is the rainfall 
rate. The dynamics of water are illustrated in Fig. [TJ see 
the left-most pair of cells (A) and the pair of cells in the 
middle (B). 



2.1.3. Plant dynamics 

We model individual plants with an IBM. Individual 
plants Ii in cell i are represented as discrete entities whose 
dynamics are given by the stochastic transition rules 



Ii — > Ei, 



3 G N{t) 



(5a) 
(5b) 
(5c) 



Here the symbols over the arrows refer to the probabilities 
per unit time, or rates, of the corresponding transitions. 
These rates may depend on the amount of available soil 
water uj (see Eqs. (O below). The first transition rule, 
Eq. (|5ap . corresponds to an individual plant in cell i 
giving birth to another plant in the same cell i at a rate 
TtiuJi). The second rule, Eq. ([5b]) . refers to the death 
of a plant in cell i. The quantity Ei on the right-hand 
side of this reaction stands for a vacancy in cell z, so that 
this transition rule indicates that an individual plant Ii 
in a cell i dies at a rate Td and leaves behind an empty 
place Ei in cell i. The third rule, Eq. (|5c|) . describes a 
process in which an individual plant in cell i gives birth 
to another plant in a neighbouring cell j at a rate T s (uij). 
This captures the dispersal of the seeds of an individual 
plant to neighbouring cells and includes the probability of 
germination of the seedling. Notice that the rate of such 
dispersion processes depends on the amount of soil water, 
ojj, in the cell in which the new plant is born (see Eq. ([B]) 
below) . 

We define the transition rates in analogy with previou s 



e.g. 



HilleRisLambers et al 



deterministic models 
IRietkerk etaD (<2002h . and IPuevo etal 
cally, we use 

r 6 (wj) = c/3(uii), 

F d = d, 
T s (uJ :i )=Kcp{uj j ) 



(120011) 

Specifi 



(6a) 
(6b) 
(6c) 



The plant birth rate Tb is proportional to the plant water 
uptake f3 (Eq. (j?])), and c is a constant that characterizes 
the conversion rate at which water uptake is turned into 
newly established plants (Eq. (JS^,)). The death rate Yd is 
assumed to be a constant d (Eq. (JBb))- The rate T s with 
which new plants are established in neighbouring cells is 
proportional to the water uptake /3, to the parameter c 
and to the constant K, representing the probability that a 
seed falls into a neighbouring cell and manages to survive 
(Eq. (Jit)). As described above the water uptake rate (3, 
and thus T s , depend on the soil water concentration in the 
colonized cell (uij). This reflects the fact that the prob- 
ability of germination of the seedling in the neighbouring 
cells is a function of the local availability of soil water. In 
our simple model we assume tha t plants can coloni ze only 
nearest neighbour cells (see e.g. iKefi et all . 2007bl ). This 



3 



can be generalized to take into account longer-range dis- 
persal kernels, which may be appropriate when the charac- 
teristic scale of seed dispersal is appreciably greate r than 



h, the lateral extension of a cell ( Puevo et al. . 20081) . It is 



worth pointing out that the transition rates do not depend 
on the whole history of the system, but only on its cur- 
rent state. This type of stochastic hybrid model is usually 
referred to as a piecewise deterministic Markov process in 
the literature (IPavisI Il984 iFaggionato et all . 120091 12010| ; 
Realpe-Gomez et all 120121 ). 



2.1.4- Coarse-graining: cell dynamics 

Given that we neglect the spatial structure within a cell, 
the relevant transition rules involve the state of cells rather 
than that of particular individual plants, so our model is of 
an effective nature; it provides a coarse grained description 
of the underlying dynamics. According to the rules for 
individual plants, Eq. ([5]), the only possible transitions 
for a cell i with m plants are rii — > rii ± 1, i.e. the birth or 
death of a plant in cell i. Hence, the possible transitions 
that can occur in cell i are ful l y spec ified by the expressions 
(see e.g. Black and McKanel ( 2012 ) for similar models) 



T b (ni + l\m; Ui) = niT b (uJi), 
Td{ni - l\ni) = n % T d , 



(7a) 
(7b) 



T' s , 3 (rii,nj + l\ni,rij; u>j) = mT s (ujj), j 6 N(i), 

(7c) 

for all i. The notation T b (ni + l\ni; Wj), for instance, stands 
for the probability per unit time (or rate) that the num- 
ber of plants ni in a cell i increases by one, given that the 
amount of soil water in i is u>i at the time the transition 
takes place. Since each plant within a cell i has the po- 
tential to give birth to a new plant in i with a rate Tb{u>i) 
independently of all other plants, the total rate for a tran- 
sition rii — > ni + 1 in cell i is given by riiY^iwi). Similar 
considerations apply for the other two processes described 
in Eq. © above. Notice that the arguments before the 
vertical bar in the above rate functions, Eqs. ([7]), refer to 
the state of the system after the transition, while those 
to the right of the vertical bar refer to the state before 
the transition, following the standard notation for condi- 
tional probabilities. This stochastic dynamics of cells is 
illustrated in Fig. [TJ see the pair of cells in the middle (B) 
and the three right-most pairs of cells (C, D, E). The cells 
at the top (C) correspond to an on-site birth event, those 
in the middle (D) correspond to a death event and those 
at the bottom (E) correspond to a birth event induced by 
a neighbouring cell. 

2.2. Deterministic approximation 
2.2.1. General considerations 

The dynamics of plants in our model are stochastic, so 
we cannot predict the exact state of the system at a future 
time. At best, it might be possible to obtain the probabil- 
ity of finding the system in a given state, specified by the 



number of plants, nj, and the soil and surface water depths 
u>i and (Tj, in all cells i. From such probability distributions 
it would then be possible to compute the expected aver- 
age behaviour. Another way to study stochastic systems 
is to perform a certain number S of independent stochas- 
tic simulations. In doing so we could, for each cell i and 
for a given time t, estimate the average density of biomass 
or the average depths of soil and surface water. We will 
denote these quantities by Pi, Wi and Oi, respectively. 

The transition rules governing the dynamics of plants 
(Eqs. (©-0) along with the dynamical laws for the soil 
and surface water (Eqs. (HJ-Q) induce a corresponding 
deterministic equation for the dynamics of these averages 
in the limit of a very large number of simulations, formally 
S — > 00. At the same time this deterministic approxima- 
tion also describes the behaviour of the stochastic model 
for very small \x. The equivalence is formally exact in 
the limit \i — > (cf. iRealpe-Gomez et al. , 2012 ) . In this 
sense fi controls the strength of the noise in the stochastic 
model: if fj, — > the noise is irrelevant and the dynam- 
ics becomes deterministic. Since pi = fj,m the number of 
plants tends to be very large (rii ~^ °o) in this limit in 
regions with a non-zero density of biomass, pi > 0. This 
indicates that a deterministic model might not be a faith- 
ful representation of a semi-arid ecos ystem with relatively 



few pl ants per cell. We also refer to Realpe- Gomez et al 



(2012J) for mo re details of a non-spatia l version of the 
model; see also Black and McKanel (|2012l) for a general re- 
view of individual-based models and the connection with 
deterministic limiting descriptions. 

The deterministic model governing the average be- 
haviour of the stochastic model turns out to be sim- 
ilar to models based on partial differential equations . 
well known in t h e liter a ture dHilleRisLambers et aL L 2001 



Rictk erk et ail , l2002t IPuevo et all 120081 iDakos et al 



20111) . The average biomass density, Pi, and the average 
depths of soil and surface water, Wi and Oi, respectively, 
follow the equations 



dP 

= F p (P i ,W i ,O i )+D p (W i ) AP { , 



dt 
dO, 



F^WuO^ + Du&Wi, 



= F a {Pi,Wi,Oi)+D a AOi. 
Here F u and F a are given by Eqs. ®-(|3]) and 

F P (P, W, O) = c/3 (W) P - dP, 



(8a) 
(8b) 
(8c) 

(9) 



where we have defined a rescaled conversion rate c = 
(1 + z K)c, with z the number of nearest neighbours of 
a cell, e.g. z — 2 in one dimension and z = 4 for two di- 
mensions. The diffusion operator is given by Eq. @ as be- 
fore, and so Eq. (I5al) above is effectively a diffusion equa- 
tion with a diffusion coefficient D p {Wi) = h 2 Kcj3 (Wi). 
The biomass diffusion coefficient thus depends on the av- 
erage amount of soil water in the cell under consideration. 



4 



This term constitutes the only difference between the de- 
terministic equations (J5J , and a discretization of the se t 
of equations introduced bv IHilleRisLambers et ail (|200lh . 
Si milarly, Eqs. dHt c oincide with the equations discussed 
in lPuevo et al.l ( 20081) in the case in which the range of the 



kernel of seed dispersal is so short that only dispersal to 
neighbouring cells matters. 

2.2.2. Homogeneous fixed points 

The deterministic equations introduced above, Eq. ©, 
have homogeneous fixed points, i.e. equilibrium states, 
such that Pi = P, Wi — W and Oi = O, for all i, and 



dP 
dt 



= 0, 



dW 
~dt 



= 0, 



dO 
dt 



= 0. 



(10) 



Such fixed points may or may not be stable under small 
spatially homogeneous perturbations. Depending on the 
choice of parameters there c an be either one o ut of two 



possible stable fixed points (jPuevo et al.l . [2008) or none 



(|Realpe-Gomez et all 120121) . The two types of stable fixed 
points correspond to a bare soil state (also referred to as 
a 'desert state' in the following) 

P d = 0, W d = -, O d = -^-, (11) 
r a W 

or to a state with non-zero vegetation 

c_R_ refer _ rffa _R_ 
a cb — a cb — d ct(P v ) 

In the region where there are no stable fixed points, a nu- 
merical integration of the corresponding non-spatial equa- 
tions shows the existence of homogeneous limit cycles 
( Realpe-Gomez et all 120121 ). 



2.2.3. Spatial patterns and resilience 

Besides the homogeneous states given by Eqs. pTp -([T2 ^) 
above, the deterministic system defined b y Eqs. flU can 



also display stationary spatial patterns ([Rietkerk et al 
2002). This is not unsurprising, given the similarity of 
Eqs. © t o those inves t igated by IHilleRisLambers et al 



( 2001 ) and Puevo et al. ( 20081) . which are known to dis- 
play spatially patterned solutions. For certain values of 
the parameters, these patterns can 'emerge' out of the ho- 
mogeneous vegetated state given by Eq. (Tj"2|) when a small 
heter ogeneous perturbation is applied dRietkerk et all 
2002 ). This results in T uring patterns (|Turintd . 1 1952 : 
Cross and Greensidel 120091) and can be studied mathemat- 
ically within a linear approximation of Eqs. (|SJ around 
the fixed point with homogeneous vegetation. Numerical 
integration of the full set of nonlinear equations has re- 
vealed, however, that spatial patterns exist in a region of 
parameters much wider than that explained by the Turing 
mechanism (Riet kerk et al ., 2002). In a parameter range 
relevant for dryland ecology (i.e. at low rainfall values) a 
relatively large region of bistability between spatial pat- 
terns and the desert state exists. In this region, and de- 
pending on the initial conditions from which the dynamics 



are started, the deterministic system can converge either 
to the desert state given by Eq. (TlTj) or to a stable vege- 
tation pattern. The basin of attraction of the state with 
vegetation patterns, i.e. the set of all initial conditions 
which lead to such a state, can be used to characterize the 
resilience of the deterministic model: the larger the size of 
the basin of attracti on of a ' desirable' state the more re- 
silien t is the system (Holli"n3. 1973 : Grimm and Calabrese , 



1 2 1 it : iMartin and Calabrese , 2011 ). Performing such a 
study quantitatively is not an easy task given that the 
space of all possible initial conditions is huge. We will 
therefore restrict our work to an analysis involving only 
one type of realistic initial conditions. This will be dis- 
cussed in more detail below. 

3. Analysis 

All results discussed in this paper were obtained from 
spatially explicit numerical simulations, mostly in one di- 
mension, with a few two-dimensional examples for illus- 
trative purposes. We used periodic boundary conditions 
throughout. In order to obtain sample simulations for the 
stocha stic mode l we us ed an ad aptation o f an a lgorithm 
due to Gillespie (Il976l) (see also Gillespie!. Il977l) . details 
are discussed in iRealpe-Gomez et al.l (|2012f ) . To obtain 
solutions of the deterministic approximation, Eq. (JSJ), we 
integrated these equations numerically using the Euler- 
forward method with a time step of At = 0.01 days. 

In all cases we used the same initial conditions for both 
the stochastic model and its deterministic approximation 
obtained as follows: (a) soil water u)i, and surface water 
Oi in all cells i were assigned the values corresponding to 
the homogeneous desert state given by Wd and Od in Eq. 
(fTTj) : (b) biomass density was assigned a value po > in a 
fraction / of the cells picked at random, and a value of zero 
in the remaining fraction (1 — /) of cells. We emphasize 
that in each run the random choice of populated cells is 
the same in both models. Notice that po, along with the 
parameter p, fixes the initial number of individual plants 
per initial populated cell no — po/ P, which has to be an 
integer. The possible choices of po and p are therefore 
constrained. 

First we analysed the differences between the determin- 
istic and the stochastic approach comparing the behaviour 
of the system in single runs. Then, in order to study 
how generic the observed outcomes were, we estimated 
the probability of extinction, P ex t, in both the stochastic 
and the deterministic models as a function of the differ- 
ent model parameters. This probability was obtained from 
simulations of the corresponding dynamics for several dif- 
ferent initial conditions, picked at random as described 
above, and counting the fraction of samples in which the 
systems became extinct. In the case of the deterministic 
system the only difference between simulation runs is the 
randomly chosen initial spatial distribution of the vege- 
tated cells. Once the initial condition is fixed, the sub- 
sequent dynamics is fully deterministic, with no further 
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random elements. Similarly, the simulation runs of the 
stochastic system also involved a random element due to 
the initial conditions, but further stochasticity then en- 
tered during the actual run due to the random processes 
determining the sequence and timing of the transitions of 
the individual-based model. 

We estimated the probability of extinction in terms of 
what we expect to be the relevant parameters of the set 
of initial conditions: (a) the initial cover /, i.e. the frac- 
tion of cells with po > in the initial conditions; (b) the 
initial amount of biomass po introduced in those cells; (c) 
the mass of a plant per unit area p; and (d) the rainfall 
rate R. We chose / and po as relevant parameters because 
they determine the amount of initial biomass, p because 
it controls the strength of the stochasticity in the model 
and R in order to investigate the models under different 
conditions of water availability. In the limit p, — > 0, the 
stochastic results are expected to show similar behaviour 
to the deterministic approximation, while for larger values 
of p we would expect stochastic effects to be more com- 
mon. 

To estimate the probabilities of extinction, 50 simula- 
tions were run for each model and for each set of parame- 
ters. Each of these simulations was run for a period of time 
of up to T = 5000 days. Simulation runs of the stochas- 
tic model were identified as becoming extinct if all cells i 
reached a state with m = within this time. In the deter- 
ministic model we applied a threshold criterion Pi < e, i.e. 
a given cell is assumed to contain no plants when its plant 
density is below a set threshold e. In order to be consistent 
with the original stochastic model we chose e = p, as we 
assumed that p was the contribution of a single plant to 
the biomass density. The introduction of a threshold was 
necessary because the deterministic approximation works 
with a genuinely continuous density, and the state with 
exactly zero biomass is reached only asymptotically at in- 
finite times. We carried out some consistency checks, and 
observe that results depend little on the exact choice of 
the value of the threshold. Findings remained essentially 
the same even for e much smaller than p. 

4. Results 

4--1- Dynamics 

We compared the dynamics of the stochastic model with 
the dynamics of the corresponding deterministic approx- 
imation. Both models were initialized with homogeneous 
initial conditions (i.e. / = I) with po = 10 g m~ 2 . For 
all cells, soil and surface water took values corresponding 
to the desert state, given by Eq. (|11[) . We compared the 
time dynamics of the total biomass density for both mod- 
els. As expected from the stability analysis of the deter- 
ministic non-spatial model, we observed that in this case 
the deterministic approximation converges asymptotically 
to the desert state (Fig. [5J red thick line). In contrast 
the stochastic model did not reach extinction for the same 



homogeneous initial conditions and parameter values. The 
stochastic model initially followed a path close to the de- 
terministic approximation and approached extinction, but 
then some regions displayed an explosive growth and the 
system finally was found in a quasi-stationary pattern. As 
expected the stochastic dynamics deviated substantially 
from its deterministic approximation especially when the 
number of plants in the system was small (see Fig. [5]). 
We will refer to this effect of escaping a path to extinc- 
tion as 'self-recovery'. From a mathematical perspective 
th is phenomenon app e ars si milar to the one investigated 



Black and McKanel ( 2011 ). For clarity the figure shows 
the results up to time t = 500 days even though we ran the 
simulations for much longer (up to time T = 5 x 10 4 days). 
We observed that the total biomass in the stochastic model 
remains essentially the same up to stochastic deviations, 
i.e. it reaches a stationary state. Still there was a nonzero 
probability that the stochastic model, after having escaped 
the initial path to extinction, e.g. after, say, t « 300 days 
in Fig. [5J reached extinction. Such events occur because 
of large stochastic deviations, and the pro bability of such 
events can be expected to be rather small ( Frev . 2010). 

The phenomenon of self-recovery in the stochastic model 
can also be seen in a comparison of the time-dependent 
spatial biomass density profile along the line of cells in 
both models, see Fig. [3] Initial conditions for the stochas- 
tic and the deterministic models were chosen as / = 1/2 
and po — 10 g m~ 2 . For all cells the soil and surface water 
densities were initialised from the values of the desert state 
given by Eq. (jlip . In Fig. [3] we compare the dynamics 
of the vegetation profile in the stochastic model with that 
of the corresponding deterministic approximation. Again, 
we can see that the deterministic approximation leads to 
extinction while the stochastic model recovers. If we look 
at the stochastic simulation for a longer time (5 x 10 4 days) 
we clearly see that vegetation persists, see Fig. 2]below the 
horizontal dashed line. We have also used the last configu- 
ration reached by the stochastic model (horizontal dashed 
line in Fig. [4]) as initial conditions for the deterministic 
model; the corresponding time dynamics is shown in Fig. 
HI above the horizontal dashed line. The patterns also re- 
main stable under the deterministic dynamics, but it was 
demographic stochasticity which allowed the system to es- 
cape extinction in the first place; running the deterministic 
dynamics alone leads to extinction. The stochastic system 
thus generated suitable initial conditions for the determin- 
istic system to converge to a spatial pattern. In Fig. 2] one 
can also observe that the patches reached by the deter- 
ministic model are fully frozen at long times while those 
generated by the stochastic model remain dynamic, they 
can split, merge or become extinct. 

Similar observations can also be made in the two- 
dimensional system, as we show in Fig. [3] and in the en- 
closed supplementary video animation. As expected, both 
the stochastic and the deterministic models display spatial 
vegetation patterns, similar to those observed previously 
with other models, and in an analogous range of rainfall 



G 



values dHilleRisLambers et alll200ll:lRietkerk et all . 12002 : 
IPuevo et alll2008h . 

Finally we notice that in many simulations the stochas- 
tic system follows a regime of 'explosive' local growth soon 
after it has escaped the path to extinction. During this 
relatively short period of time some cells can get to carry 
a relatively large number of individuals with a maximum 
of about 120 plants observed in simulations. Afterwards 
plants spread out in space to a certain extent, and even- 
tually the system appears to reach a stochastic stationary 
state (see supplementary video). 

4-2. Probability of extinction 

Next we will present the results of a more systematic in- 
vestigation of the effect of 'self-recovery'. Figure [5] shows 
the probability of extinction P cxt as a function of /, for 
two different values of po (10 g m~ 2 and 50 g m -2 ). We 
observe that the stochastic model almost never reached ex- 
tinction, while the deterministic approximation did so in 
almost all samples with / > 0.4. These observations may 
appear counter-intuitive at first sight: in the deterministic 
model we find that the more plants are in the system ini- 
tially, i.e. the larger the initial cover /, the more likely it is 
that the system reaches extinction. However, this effect is 
already seen in the deterministic non-spatial model which 
predicts the extinction of homogeneous initial configura- 
tions (i.e. / = 1). The spatial model instead predicts sta- 
ble pattern s ( HillcR isLambers et al. , 2001 ; Rietkerk et al 



2002- Puevo et al.l . l2008l) . Even more surprisingly perhaps, 



is that the larger the initial amount of biomass in each 
populated cell, po , the larger the probability of extinction 
tends to be. From the simulations, it seems that if the 
amount of biomass in a cell is too large, a fast spatial 
spread of plants is triggered, which favours a more homo- 
geneous distribution of biomass. This amounts to an effec- 
tive increase of the initial cover, /, and thus the previous 
effect takes over. 

We study the variations of the probability of extinction 
Pext as a function of ji for three different values of /: 1/8, 
1/2, and 7/8 in Fig. [7] We can observe that P ex t ~ 
for almost all values of n except for // = 0.1 g m~ 2 and 
/i > 6 g m~ 2 , which correspond to very small and large 
plant biomass, respectively. 

Next, we study the probability of extinction P oxt as a 
function of the rainfall rate R, for a range corresponding to 
typical dryland values, using three different values for /1: 1 
g m -2 , 5 g m~ 2 and 10 g m~ 2 (Fig. [3]). In the stochastic 
model, the effect of 'self-recovery', or escaping the path 
to extinction, is more appreciable the larger the amount 
of rainfall. The corresponding probability of extinction 
for the deterministic system is equal to one in almost the 
whole regime investigated (not shown). 

Under certain conditions, the opposite effect can hap- 
pen, with demographic noise promoting extinction in cases 
in which the deterministic system survives instead. How- 
ever, this needs a range of parameters which we expect not 
to be relevant for real-world situations. In particular, we 



have observed such behaviour only for large values of the 
parameter /i and for low initial cover. More specifically, 
we have observed that this can happen if / < 10% and 
\i > 10 g m~ 2 , for R > 0.5 mm d _1 , or /i > 5 g m~ 2 , for 
0-4 < R < 0.5 mm d -1 . 

5. Discussion 

Demographic stochasticity is known to have impor- 
tant consequences that are usually neglected when mod- 
cllin g ecosystems in te r ms o f continuous densities (see 
e.g. Black and McKanel 2012 ). In this work, we repre- 
sented semi-arid ecosystems with two parallel models: a 
hybrid stochastic model, including individual-based vege- 
tation dynamics and deterministic water dynamics, and a 
deterministic model, governing the behaviour of the aver- 
age stochastic variables. The models represented explicitly 
the water-vegetation infiltration feedback. Thus, for the 
parameter values we studied, both models contained solu- 
tions corresponding to stable vegetation patterns, analo- 
gously to what was observed p reviously with other m odels 



of this type of ecosystem (e.g. iRietkerk et all 12002ft 



Our analysis showed that the resilience of the vegetation 
patterns changed when taking into account the discrete 
nature of plants and the intrinsic stochasticity of their 
behaviour. Including such stochasticity, the modelled 
ecosystems did not turn into deserts in a wide range of 
cases in which the deterministic representation would 
predict the ecosystem to become extinct. This is an 
important observation, given that semi-arid ecosystems 
are characterized by a rather scarce number of plants, 
scattered across regions of empty land. In a finite pop- 
ulation, the number of individuals varies because of the 
intrinsic stochastic behaviour of the individuals. Usually, 
the fewer the indi viduals, the stronger the effects of the de 
mogr a phic noise dRamaswamv et aL L 2012 1_ Rogers et al 



2012; 



2011 



Biancalani et al 



Biancalani et al 



2011 



2010; 



Butler and Goldenfeld 



Butl er and G oldcnfcld 



2009t iMcKane and Newmanl . l2005h~ Intuitively, this is 
expected to promote extinction when the number of 
plants is small. It is therefore remarkable that we ob- 
served a relevant regime of parameters in which including 
the stochastic individual nature of the plants actually 
increased the likelihood of vegetation pattern emergence, 
and of escaping the desert state. The opposite might also 
occur, i.e. a larger probability of becoming extinct in 
the stochastic model, especially when considering large 
individual plants, covering initially only a small fraction 
of soil. Under these conditions, the probability of a large 
stochastic deviation that brings the system suddenly to 
extinction, could be non- negligible. We expect, however, 
that our approach is better justified for intermediate 
values of [A, i.e. the biomass of a plant (per unit area), 
where the results were rather robust. When the individual 
plant biomass was relatively high (/x > 8 g m~ 2 ), or very 
low (/j, — > 0), the outcomes of the stochastic model 
tended to have higher probability of being desert than for 
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intermediate /x values. Demographic noise seemed to be 
more important in an intermediate range of single-plant 
biomass, which correspo nds to a realis tic range for herbs 
or grass biomass (see e.g. iPetersl . l2002h . 

In order to study the resilience of the vegetation pat- 
terns, we addressed the issue of how and how easily 
they were reached in time, i.e. for which type of initial 
conditions the system evolved towards the pattern (e.g. 
lEppingal 



2009). This gave an indication of how "attrac- 



tive" the patterns were. In the case of the deterministic 
approximation, this was equivalent to studying the basin of 
attraction of the patterned stable states. This notion, how- 
ever, is not applicable to stochastic systems, where there 
is no notion of equilibrium point, and no deterministic dy- 
namics uniquely leading from a certain initial condition to 
a final state. For this reason, we focused on the probability 
of extinction as a function of the initial conditions and for 
different parameter values. For a large set of initial con- 
ditions and parameters, we observed that the stochastic 
model almost never led the vegetation to extinction, while 
the vegetation in the deterministic model almost always 
became extinct. When varying rainfall in a realistic range 
for drylands, the effect of 'self-recovery' in the stochastic 
model was promoted for the highest rainfall values, where 
the contrast with the deterministic case was sharper. 

The differences between the two models were especially 
important for vegetation initially covering more than 30- 
40%, where the deterministic case would practically always 
become extinct. This appeared counter- intuitive, in the 
sense that a larger initial number of plants enhanced ex- 
tinction. The more homogeneously distributed the plants 
were, the more difficult it was for any particular vegetation 
patch to actively increase its own soil water. The water- 
vegetation feedback needed contrasts between vegetated 
and bare patches to take place effectively. For the same 
reasons, in this regime the homogeneous vegetated state 
was not stable, and full homogeneous vegetation cover 
would quickly evolve into a desert state. On the other 
hand, for very low initial vegetation cover (less than 20- 
30%), the results of the deterministic and the stochastic 
approaches were very similar. We must underline that 
these low fractions of vegetation cover are quite realistic 
in the most arid ecosystems. Therefore, the deterministic 
model studied in this paper seemed to be more accurate at 
the driest end of the range considered, where initial vege- 
tation cover and rainfall were low. One may wonder how 
it is that for large initial cover (e.g. larger than 40%) the 
deterministic system has a large probability of extinction, 
while for small cover (e.g. smaller than 20%) it almost al- 
ways survives; after all the dynamics is continuous and to 
reach extinction the system first needs to go through states 
of small cover. However, the spatial structure of the inter- 
mediate states thus reached does not necessarily coincide 
with that of a state with homogeneous distribution of soil 
and surface water and a random distribution of constant 
biomass po, as were used here as initial conditions. 

The deterministic model we investigated in this work 



describes exactly the typical behaviour of the stochastic 
system when the individual plant biomass was negligible 
(/i — > 0). In this case we could have a very large num- 
ber of plants per cell. This deterministic model was a 
good approx i matio n to the kernel model investigated by 



Puevo et al.l (|2008h in the case that the range of the ker- 



nel of seed dispersal in the latter was of the order of one 
or two cells (~ 2-4 m). Under these conditions, the kernel 
could be approximated with an effective diffusive term and 
with a diffusion coefficient depending on the state of the 
soil water. If, additionally, we could neglect the spatial 
variation of the soil wat er, the two models became sim ilar 
to the model studied by iHilleRisLambers et al. ( 2001 ). In 
this sense, we could say that the stochastic model intro- 
duced here was close to these well-studied models and, in 
particular, it was not unexpected to observe spatial pat- 
terns in a similar range of parameters. What could change 
drastically is the transitory dynamics while reaching a sta- 
tionary state, as we showed with the effect of self-recovery. 

However, we should admit that it is not completely clear 
whether plants are better described as unitary individuals 
to be count ed, as we did, or rather as a continuous density 



of biomass (|Crawlevl . Eflflflh . Plants are extremely plastic 



and, unlike animals, can partially die and recover later on. 
In principle, a more complete approach would have to deal 
with each single plant individually, with all its attributes 
and ongoing processes, as has been done for instanc e in the 
study of forests, see e.g. IPerrv and Enright (l200fih . How- 
ever, this would make the problem far more complex from 
a computational point of view, and it might become diffi- 
cult to gain insight. Indeed, as has been discovered in the 
investigations on forests, a far simpler analytical approxi- 
mation, called the perfect plastic approximation, may cap- 
ture most of th e relevant details of the dynamics observed 
in simulations (jStrigul et all [2008). This case, however, 



corresponds to a regime of high vegetation density, where 
fluctuations of the average behaviour are expected to be 
irrelevant. This raises the question of which would be the 
'best' way of modelling a plant in the study of semi-arid 
ecosystems in the sense of, paraphrasing Einstein, 'keeping 
it simple, but not too simple'. 

In our stochastic modelling approach, we did not con- 
sider the heterogeneity in plant size. In particular, we did 
not include the different plant life stages. In a sense, plants 
were instantaneously created as adult individuals. This 
might appear as an unrealistic feature that might promote 
the effect of 'self-recovery' we discussed, because a plant 
could start increasing the availability of water locally as 
soon as it was created, and therefore instantly start pro- 
moting its own survival. However, this feature would also 
lead to an overestimation of water uptake. Since in our 
model plants died suddenly, in detriment of self-recovery, 
mortality was also over-estimated. Furthermore, in the 
stochastic model we investigated, a plant could produce 
another plant only in the neighbouring cells, thus limiting 
the impact the new birth had on the state of the ecosys- 
tem. 
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A next step in the complexity of the modelling could 
be to introduce two type of individuals: seedlings and 
established plants. In this way we would be able to 
take into account, for instance, the high asymmetry in 
mortality between these two. The question is then how 
robust are the results we have discussed in this work 
under this more realistic scenario. At first sight, one could 
expect that the stochastic model becomes closer to the 
deterministic approximation, but experience in this field 
of research has shown that counter-intuitive effects are 



not u n common ( Ramaswa mv et all 2012: Roger s et al 



2012; 



2011 



Biancalani et al. 



Biancalani et all 



2009; Mc Kane and Newma n. 2005). 




Butler and Goldcnfcld 
Butler and Goldenfeld 



The model we presented did not include environ- 
mental heterogeneity and stochasticity, or topography. 
We also discarded rainfall intcrmittence. Vegetation 
in arid areas is well adapted to the occasional occur- 
rence of rainfa ll (Baude na et all 2007; Kl etter et al. , 2009: 



D'Odorico et all 120071) ." although the effect of rainfall 



intermittency may not be to o relevant when spatial 
feedb acks are represented (see Baudena and ProvenzaleL 
2008). Besides water, we did not consider any other 



limiting resource, such as nutrient or light limitations. 
Moreover, we did not include another water-vegetation 
feedback mechanism, which is known to play a role in 
drylands, namely the effect of root length. Plant wa- 
ter availability increases with the root extent, which in 
turn increases with the biomass itself, thus fa vouring plant 
growth and generating vegetation patterns (iGilad et all 
2004tlLefever. R and Leieunei . ll997tlBarbier et alll2008l ). 



Another relevant issue is how to validate our findings 
with observations. For example, we could analyse patch 
dynamics to see whether it displays the 'wiggly' behaviour 
observed in the model results. This nevertheless would re- 
quire time series of spatial patterns, with an appropriate 
spatial and temporal resolution, which may not be easily 
obtainable. In principle, it could even be possible to com- 
pute the statistics of this patch dynamics, in order to have 
quantitative predictions. 

Despite its limitations, this investigation indicated that, 
in certain regimes, including demographic noise could lead 
to a larger estimate of the resilience of semi-arid ecosys- 
tems. Our model results suggested that demographic noise 
may be more important in the less arid ecosystems, with 
higher rainfall and vegetation cover. Since changes in rain- 
fall regimes are expected, for example as a consequence of 
climate change, it may be necessary to take into account 
individual-based dynamics to evaluate the resilience and 
resistance of these ecosystems to such forcing. 
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STOCHASTIC BIRTH OR 
DEATH OF PLANTS AT t + r 



DETERMINISTIC EVOLUTION 
OF WATER BETWEEN t AND t + r 

Time t Time t + r 

Eqs. ©-gl) 



Hi 


rij 




Hi 




| n,T d 


Hi - I 








Ui, (Ti 


Uj, (Tj 




LOi, Oi 



A 



B 



rii + 1 

Ui, Oi 



D 



LOi, Oi 



n i + 1 

Uj, Oj 



niT B (u)j) 



Figure 1: (Colour online) Illustration of the stochastic hybrid model for semi-arid ecosystems. For clarity we show only two neighbouring cells 
i and j but the same applies to all the other cells. Suppose that, at time t, the system is in a state where the number of plants and the depths 
of soil and surface water in cell i are given by n;, and a^, respectively (A). The analogous quantities in cell j are n^uj'yjj'y Suppose, 
furthermore, that the next birth or death of a plant takes place at time t + r; these events happen at random, and so r is itself a stochastic 
variable. Since there are no transition events between t and t + r, soil and surface water in all cells evolve deterministically according to Eqs. 

in this time interval. Suppose now that their new state at t + r is given by uij and at, for all cells i (B). At t + r a stochastic transition 
happens, there are three possible types of transitions: a plant at a cell i gives birth to a new plant in the same cell (C), it dies (D) or it gives 
birth to a new plant in a neighbouring cell j (E). These transitions happen with rates n^r^u^), n^T^, and niT s (u)j), respectively (see Eqs. 
(0 and 10). Immediately after t + r, and until the next birth or death of a plant takes place, soil and surface water in all cells again evolve 
deterministically as before. 
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Figure 2: (Colour online) Comparison of the dynamics of the to- 
tal biomass density in the stochastic model (continuous green line) 
with that in the corresponding deterministic approximation (thick 
red line). While the deterministic approximation leads to extinc- 
tion, the full stochastic model recovers. Simulations were run on a 
line with 128 cells, periodic boundary conditions and the same ran- 
dom initial conditions in all cells: soil and surface water depths were 
given by the values in the desert state, Eq. {TTJ, while the initial 
biomass was po = 10 g m — 2 (/ = 1). Other key parameter values 
were: /i = 1 g m — 2 , R = 0.6 mm d — 1 . See Table^for the remaining 
parameter values. 
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Figure 3: (Colour online) Comparison of the dynamics of the vegeta- 
tion profile in the stochastic model (left) with that in the correspond- 
ing deterministic approximation (right). While the deterministic ap- 
proximation leads to extinction, the full stochastic model recovers. 
Simulations were run on a line with 128 cells, periodic boundary con- 
ditions and the same initial conditions for both the stochastic model 
and its deterministic approximation: soil and surface water depths in 
all cells were given by the values in the desert state, Eq. 1111 . while 
biomass was po = g m~ 2 in half of the cells, chosen randomly, and 
po = 10 g m -2 in the other half (/ = 1/2). Key parameter values 
were: p = 1 g m — 2 , R = 0.6 mm d — . See Table^for the remaining 
parameter values. 



Figure 4: (Colour online) Dynamics of the vegetation profile pre- 
sented on the left side of Fig. \3\ (stochastic model) for a longer 
period of time. Below the horizontal dashed line, i.e. from t = 
days to t = 5 X 10 4 days, we show the dynamical behaviour of the 
stochastic model, and observe that the system indeed escapes the 
path to extinction and finally reaches a stationary pattern. Vege- 
tation patches, however, appear to follow dynamics on their own: 
they can split, merge, and become extinct. In order to compare 
this with the dynamics of the deterministic model, we use the con- 
figuration reached by the stochastic dynamics at t = 5 X 10 4 days 
(horizontal dashed line) as an initial condition for the deterministic 
model. The outcome of the corresponding deterministic dynamics 
are shown above the horizontal dashed line, i.e. from t = 5 X 10 4 
days to t = 7 X 10 4 days. In other words, at time t = 5 X 10 4 days we 
switch the dynamics from the stochastic to the deterministic model. 
Clearly the patterns remain stable under the deterministic dynam- 
ics. Simulations were run on a line with 128 cells, periodic boundary 
conditions and the same initial conditions for both the stochastic 
model and its deterministic approximation: soil and surface water 
depths in all cells were given by the values in the desert state, Eq. 
mil l, while biomass was po = g m -2 in half of the cells, chosen 
randomly, and po = 10 g m~ 2 in the other half (/ = 1/2). Key 
parameter values were: fi = 1 g m~ 2 , R = 0.6 mm d — 1 . See Tablc[TJ 
for the remaining parameter values. 
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Figure 6: (Colour online) Probability of extinction in the stochastic 
model (blue triangles and magenta rhombi) and in the corresponding 
approximation to a deterministic model (red squares and green cir- 
cles). We can observe a very strong contrast between the two: while 
the deterministic model almost always becomes extinct for cover val- 
ues / > 0.4, the stochastic system almost never does, for any value 
of /. For the stochastic case, the probabilities were estimated from 
50 numerical simulations for each point. The horizontal axis indi- 
cates the initial fraction / of populated cells, i.e. cells with initial 
biomass po > 0. Simulations were run on a line with 128 cells and 
with periodic boundary conditions. Both the stochastic model and 
its deterministic approximation were started with the same initial 
conditions: in all cells soil and surface water depths were given by 
the values in the desert state, Eq. 11U . while biomass was po = 10 
g m — 2 (green circles and magenta rhombi), and po = 50 g m — 2 (red 
squares and blue triangles) in a fraction / of randomly chosen cells 
and zero in the remaining cells. Key parameter values were: fi = 1 
g m — 2 , R = 0.6 mm d — 1 . See Table^for the remaining parameter 
values. 
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Figure 5: (Colour online) Comparison of the dynamics of the vegeta- 
tion profile in the stochastic model in two dimensions (left) with that 
in the corresponding deterministic approximation (right); the axes 
correspond to the two spatial dimensions. See also the supplementary 
video which shows the full dynamics of these same vegetation pro- 
files. While the deterministic approximation leads to extinction, the 
full stochastic model recovers. The simulations were run on a square 
grid of 64 X 64 cells and with periodic boundary conditions. Both the 
stochastic model and its deterministic approximation were started 
with the same initial conditions: soil and surface water depths in all 
cells were given by the values in the desert state, Eq. 1111 . while 
biomass was po = g m~ 2 in half of the cells, chosen randomly, and 
po = 10 g m -2 in the other half (/ = 1/2). Key parameter values 
were: p = 1 g m — 2 , R = 0.6 mm d — . See Table[T]for the remaining 
parameter values. 
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Figure 7: (Colour online) Probability of extinction in the stochastic 
model as a function of the parameter p, i.e. the average mass of a 
plant divided by the area of a cell. The probabilities were estimated 
from 50 numerical simulations for each point. Simulations were run 
on a line with 128 cells with periodic boundary conditions and initial 
conditions chosen as follows: in all cells initial soil and surface water 
depths were given by the values in the desert state, Eq. dill) , while 
biomass was po ~ 10 g m — 2 in a fraction / of randomly chosen cells 
and zero in the remaining cells. The three curves correspond to three 
different values of /: 1/8 (red squares), 1/2 (green circles) and 7/8 
(blue triangles). See TablelTlfor the remaining parameter values. For 
the same parameter values and initial conditions the probability of 
extinction in the deterministic model is essentially one throughout 
the whole regime investigated (not shown). 
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Figure 8: (Colour online) Probability of extinction in the stochastic 
model as a function of the rainfall rate R. Extinction probabilities 
were estimated from 50 numerical simulations for each point. Sim- 
ulations were run on a line with 128 cells with periodic boundary 
conditions. In all cells the initial depths of soil and surface water 
were given by the values in the desert state, Eq. 111! , while the 
initial biomass was po = 10 g m — 2 in half of the cells (/ = 1/2) 
and zero in the remaining half. The parameter /i took three different 
values for each of the three curves shown: 1 g m~ 2 (red squares), 5 g 
m -2 (green circles) and 10 g m — 2 (blue triangles). Note that for the 
same parameters and initial conditions the probability of extinction 
of the deterministic model is essentially one over the whole regime 
investigated (not shown) . See Table ^ for the remaining parameter 
values. 
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Parameter 


Description 


Units 


Value 


a 


maximum infiltration rate 


d- 1 


0.2 


b 


maximum specific water uptake 


mm g m d 


0.05 


c 


conversion of water uptake to plants 


g mm -1 m~ 2 


10 


d 


plant mortality rate 


d- 1 


0.25 


r 


water loss due to drainage and evaporation 


d- 1 


0.2 


h 


length of the side of a cell 


m 


2 




half-saturation constant of water uptake 


mm 


5 


k 2 


half-saturation constant of water infiltration 


g m _1 


5 




mean contribution of a plant to biomass density 


g m 


0.1-10 




diffusion coefficient for soil water 


m 2 d- 1 


0.1 


D a 


diffusion coefficient for surface water 


m 2 d- 1 


100 


K 


probability of a seed moving to a neighbouring cell 




0.02 


L 


number of cells 




128, 64 x 64 


R 


rainfall rate 


mm d _1 


0.35-0.60 


Wo 


water infiltration rate in the absence of plants 




0.1 



Tabic 1: Parameter values for the models studied in this paper. 
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